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Abstract.  This  paper  considers  the  comparative  eval- 

uation of  algorithms  for  mathematical  programming  problems. 

It  is  concerned  with  the  measurement  of  computational  speed 
and  examines  critically  the  concept  of  equivalent  number  of 
function  evaluations  Ng.  Does  this  quantity  constitute  a 
fair  way  of  comparing  different  algorithms? 

The  answer  to  the  above  question  depends  strongly  on 
whether  or  not  analytical  expressions  for  the  components  of 
the  gradient  and  the  elements  of  the  Hessian  matrix  are 
available.  It  also  depends  on  the  relative  importance  of 
the  computational  effort  associated  with  algorithmic  opera- 
tions vis-a-vis  the  computational  effort  associated  with 
function  evaluations. 

Both  theoretical  considerations  and  extensive  numerical 
examples  carried  out  in  conjunction  with  the  Fletcher-Reeves 
algorithm,  the  Davidon-Fletcher-Powell  algorithm,  and  the 
quasilinearization  algorithm  suggest  the  following:  the 

concept,  while  accurate  in  some  cases,  has  drawbacks  in 
other  cases;  indeed,  it  might  lead  to  a distorted  view  of 
the  relative  importance  of  an  algorithm  with  respect  to 
another. 

The  above  distortion  can  be  corrected  through  the  in- 
troduction of  a more  general  parameter  . This  generalized 
parameter  is  constructed  so  as  to  reflect  accurately  the  com- 
putational effort  associated  with  function  evaluations  and 
algorithmic  operations. 
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From  the  analyses  performed  and  the  results  obtained,  it 
is  inferred  that,  due  to  the  weaknesses  of  the  Ng  concept, 
the  use  of  the  Ng  concept  is  advisable.  In  effect,  this  is 
the  same  as  stating  that,  in  spite  of  its  obvious  shortcomings, 
the  direct  measurement  of  the  CPU  time  is  still  the  more 
reliable  way  to  compare  different  minimization  algorithms. 

Key  Words.  Numerical  analysis,  numerical  methods,  computing 
methods,  computing  techniques,  complexity  of  computation, 
philosophy  of  computation,  comparison  of  algorithms,  computa- 
tional speed,  measurement  of  computational  speed,  number  of 
function  evaluations,  equivalent  number  of  function  evalua- 
tions, unconstrained  minimization,  mathematical  programming. 
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1.  Introduction 

Over  the  past  two  decades,  a large  number  of  mathematical 
programming  problems  have  been  studied  both  analytically  and 
numerically.  Generally  speaking,  these  problems  belong  to 
four  principal  categories:  (i)  unconstrained  minimization 

problems,  (ii)  constrained  minimization  problems  involving 
equality  constraints,  (iii)  constrained  minimization  problems 
involving  inequality  constraints,  and  (iv)  constrained  mini- 
mization problems  involving  both  equality  and  inequality  con- 
straints. 

For  each  category  of  problems,  three  types  of  methods 
have  been  developed,  more  specifically:  (a)  zeroth-order 

methods,  (b)  first-order  methods,  and  (c)  second-order 
methods.  Methods  of  type  (a)  utilize  only  the  functions 
under  consideration  and  avoid  the  computation  of  derivatives. 
Methods  of  type  (b)  utilize  the  functions  under  consideration 
and  their  first  derivatives.  And  methods  of  type  (c)  utilize 
the  functions  under  consideration  together  with  their  first 
and  second  derivatives. 

For  each  category  of  problems  and  each  method,  several 
classes  of  algorithms  have  been  developed.  As  an  example, 
with  reference  to  the  category  of  unconstrained  minimization 
problems  and  first-order  methods,  the  following  classes  of 
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algorithms  are  available  today:  ordinary-gradient  algorithm, 

conjugate-gradient  algorithm,  variable-metric  algorithm, 
memory-gradient  algorithm,  and  supermemory-gradient  algorithm. 

It  is  clear  that  a bewildering  combination  of  problems, 
methods,  and  algorithms  exists  and  that  the  proliferation  of 
these  algorithms  is  bound  to  cause  some  confusion  in  the 
user,  that  is,  the  engineer,  the  chemist,  or  the  economist 
who  must  solve  problems  of  the  real  world.  Faced  with  a 
given  technical  problem,  the  user  would  like  to  know  an 
answer  to  the  following  question:  what  kind  of  algorithm 

should  be  selected  to  solve  the  problem  under  consideration? 

Unfortunately,  no  clear-cut  answer  can  be  given  to  the 
above  question.  Nevertheless,  the  identification  of  poten- 
tially successful  algorithms  can  be  facilitated  if  the  de- 
veloper of  an  algorithm  supplies  sufficient  information 
about  the  following  items:  (A)  algorithm  robustness  or 

convergence  range;  (B)  convergence  rate;  (C)  computational 
speed;  (D)  memory  requirements;  and  (E)  programming  com- 
plexity. Only  for  simply-structured  problems  (for  instance, 
linear-quadratic  problems) , the  above  information  can  be 
predicted  theoretically.  For  more  general  problems,  the 
help  of  computer  experimentation  is  nearly  indispensable. 
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In  this  study,  we  are  concerned  with  Item  (C) , the 
measurement  of  computational  speed.  In  particular,  we 
examine  critically  the  concept  of  equivalent  number  of 
function  evaluations  Ng  and  inquire  whether  this  quantity 
constitutes  a fair  way  of  comparing  different  minimization 
algorithms.  Then,  we  introduce  a more  general  parameter 
Ng  , which  is  constructed  so  as  to  reflect  accurately  the 
computational  effort  associated  with  function  evaluations 
and  algorithmic  operations.  Next,  we  examine  the  Ng  con- 
cept vis-a-vis  the  Ng  concept  through  several  numerical 
examples.  These  examples  include  some  widely  used  test 
functions  (Rosenbrock,  Wood,  Powell,  and  Miele  functions 
plus  generalized  Rosenbrock  function)  and  some  widely  used 
minimization  algorithms  (Fletcher-Reeves  algorithm,  Davidon- 
Fletcher-Powell  algorithm,  and  quasilinearization  algorithm) . 

Unconstrained  Minimization.  For  the  sake  of  simplicity, 
we  consider  in  the  following  sections  only  one  category  of 
problems,  namely,  unconstrained  minimization  problems.  More 
specifically,  we  consider  the  function 

f = f(x)  , (1) 

where  f is  a scalar  and  x is  an  n-vector  whose  components 
are  unconstrained.  We  denote  by  g(x)  the  gradient  and  by 
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H (x)  the  Hessian.  We  observe  that  the  gradient  vector  has 

2 

n components  and  that  the  Hessian  matrix  has  n elements. 

2 

Of  these  n elements,  only 


m = n(n+l)/2 


(2) 


need  to  be  calculated,  owing  to  the  symmetry  of  the  Hessian 
matrix. 


< 

L 
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2 . Measurement  of  Computational  Speed 

The  methods  employed  for  measuring  the  computational 
speed  of  different  minimization  algorithms  can  be  grouped 
into  two  classes:  (i)  direct  measurement  and  (ii)  indirect 

measurement. 

Direct  Measurement.  The  most  direct  way  for  evaluating 
the  computational  speed  of  an  algorithm  is  to  measure  the 
so  called  CPU  time  (the  symbol  CPU  stands  for  central  pro- 
cessing unit) . The  main  advantage  of  this  quantity  is  that 
it  includes  both  function  evaluation  time  and  algorithmic 
time.  The  main  disadvantage  is  that  the  CPU  time  is  machine 
dependent  as  well  as  operator  (programmer)  dependent.  The 
above  difficulties  can  be  removed  to  some  degree  if  the  com- 
parison of  different  algorithms  is  done  on  a single  computer, 
with  the  same  programming  language,  with  the  same  compiler, 
with  the  same  subroutines,  under  similar  workload  conditions 
of  the  computer,  and  by  the  same  programmer.  In  other  words, 
it  is  essential  that  the  same  experimental  conditions  be  kept 
for  all  of  the  algorithms  being  investigated.  For  an  ex- 
ample of  comparative  experiments  done  under  these  conditions, 
see  Refs.  1-2. 

Normalized  Time.  In  an  attempt  to  make  the  direct  measure- 
ment of  the  CPU  time  independent  of  the  particular  computer, 
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Colville  introduced  in  Ref.  3 the  concept  of  normalized  time: 

Tn  = T/Ts  • (3) 

Here,  T is  the  CPU  time  required  to  solve  a particular  test 
problem  with  the  algorithm  under  consideration,  and  Tg  is 
the  CPU  time  required  to  execute  a so-called  standard  pro- 
gram, devised  by  Colville.  This  standard  program  consists 
of  inverting  a 40x40  matrix  ten  times.  Ideally,  this  para- 
meter should  be  machine  independent.  In  practice,  it  still 
depends  on  the  subroutines  and  the  compiler  used  (Ref.  4). 

Indirect  Measurement.  There  exist  three  major  ways  for 
evaluating  the  computational  speed  of  an  algorithm  indirectly: 
(a)  number  of  iterations;  (b)  number  of  function,  gradient, 
and  Hessian  evaluations;  and  (c)  equivalent  number  of  function 
evaluations.  All  of  these  quantities  are  machine  independent, 
operator  independent,  and  simple  to  compute.  However,  their 
use  implies  the  drawbacks  discussed  below. 

Number  of  Iterations.  Generally  speaking,  the  time  per 
iteration  varies  from  one  algorithm  to  another.  Therefore, 
one  cannot  employ  the  number  of  iterations  N as  an  indicator 
of  computational  speed,  unless  one  can  be  reasonably  sure 
that  all  of  the  algorithms  being  compared  require  approxi- 
mately the  same  workload  per  iteration.  For  an  example  where 
this  situation  arises,  see  Ref.  5. 
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Number  of  Function,  Gradient,  and  Hessian  Evaluations.  If 
one  can  be  reasonably  sure  that  the  algorithmic  time  is 
negligible  by  comparison  with  the  function  evaluation  time, 
one  can  use  the  triplet  composed  of  number  of  function  evalua- 
tions (Nq) , number  of  gradient  evaluations  (N^),  and  number 
of  Hessian  evaluations  (N2)  as  a collective  indicator  of  computa 
tional  speed.  The  trouble  is  that  the  resulting  indication  is 
unclear,  unless  one  is  willing  to  attribute  relative  weights 
to  function,  gradient,  and  Hessian  evaluations. 

Equivalent  Number  of  Function  Evaluations.  Let  the  rela- 
tive weights  (l,n,m)  be  attributed  to  the  elements  of  the 
triplet  (Ng,  N^,  N2 ) . With  this  understanding,  one  can  form 
the  following  linear  combination: 


N = N_  + nN,  + mN„  , 

e 0 12 


(4) 


which  is  called  the  equivalent  number  of  function  evaluations. 
This  parameter  can  be  used  as  an  indicator  of  computational  speed 


providing  the  algorithmic  time  is  negligible  by  comparison 
with  the  function  evaluation  time  and  providing  the  weights 
(l,n,m)  measure  correctly  the  relative  importance  of  function 
evaluation,  gradient  evaluation,  and  Hessian  evaluation. 
However,  this  depends  on  whether  analytical  expressions  for 
the  components  of  the  gradient  and  the  elements  of  the  Hessian 
matrix  are  available  or  not. 


1 
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3 . Standard  Definition  of  Equivalent  Number  of  Function 
Evaluations 

Assume  that  some  particular  algorithm  is  employed  in  order 
to  obtain  the  minimum  of  the  function  (1)  on  a digital  computer. 
The  total  CPU  time  T can  be  written  as 

T = T + T • (5) 

a e ' ' 

Here,  T is  the  algorithmic  time,  namely,  the  CPU  time  required 
to  perform  the  arithmetic  operations  intrinsic  to  the  al- 
gorithm being  employed.  And  Tg  is  the  function  evaluation 
time,  namely,  the  CPU  time  required  to  evaluate  the  function, 
the  gradient,  and  the  Hessian.  Therefore,  Tg  can  be  written  as 


T 

e 


+ T,  + T„ 


(6) 


Here,  T^  denotes  the  CPU  time  associated  with  function 
evaluations,  T^  denotes  the  CPU  time  associated  with  gradient 
evaluations,  and  T^  denotes  the  CPU  time  associated  with 
Hessian  evaluations. 

Let  Tq,  t , t ^ denote  the  basic  times  required  to  compute  one 
function,  one  gradient,  and  one  Hessian,  respectively.  Observe 
that  the  following  relations  hold: 
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T0  T0N0  ' T1  T1N1  ' T2  T2N2  ’ 


As  a consequence,  Eq.  (5)  can  be  rewritten  as 


(7) 


T 


T + 
a 


Vo  + 


T1N1  + T2N2 


(8) 


Next,  let  the  following  assumptions  be  employed: 

(Al)  From  the  point  of  view  of  the  CPU  time,  one  gradient 
evaluation  is  equivalent  to  n function  evaluations. 

(A2)  From  the  point  of  view  of  the  CPU  time,  one  Hessian 
evaluation  is  equivalent  to  m function  evaluations. 

(A3)  The  algorithmic  time  is  negligible  by  comparison 
with  the  function  evaluation  time. 

In  equation  form,  the  above  assumption  can  be  rewritten 
as  follows: 


(Al)  = nig  / (9-1) 

(A2 ) x 2 = nix-  f (9-2) 

(A3)  T <<  T . (9-3) 

cl  6 


As  a consequence, Eq.  (8)  can  be  rewritten  as 
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T = t.N 
0 e 


where 


Ne  = Nq  + nNx  + mN2  . 


This  is  the  standard  definition  of  equivalent  number  of 
function  evaluations  for  second-order  methods,  which  reduces  to 


Ne  = Nq  + nU1 


for  first-order  methods  (N2  = 0)  , and  to 


N = Nn 
e 0 


for  zeroth-order  methods  (N^  = n2  = 0) 


_ 


* 


T" 
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4 . Cases  Where  the  Assumptions  Are  Satisfied 

In  this  section,  we  present  some  cases  where  the 
equivalent  number  of  function  evaluations  (11)  can  be  regarded 
to  be  a correct  indicator  of  computational  speed. 

Example  4.1.  Suppose  that  the  function  (1)  is  such  that 
analytical  expressions  for  the  components  of  the  gradient  and 
the  elements  of  the  Hessian  matrix  are  not  available.  Con- 
sequently, some  numerical  approximation  scheme  to  the  gradient 
and  the  Hessian  is  necessary.  For  example,  suppose  that  a 
forward  difference  scheme  is  employed.  Denote  by  e some 
small  number,  and  denote  by  u^  a unit  vector  in  the  x^-  direction. 
Then,  the  following  relations  hold: 

f(x  + eui)  - f ( x ) = eu^g(x)  = egi  (x)  , (14) 

where  i = l,2,...,n,  and 

f (x+eu^+eu j ) -f (x+eu^ ) -f (x+eu ^ ) +f (x)  = c^u^H(x)u^  = e^h^  (x)  , (15) 
where  i = 1 , 2 , . . . , n and  j = i ,....,  n 


If  f(x)  is  known,  Eq.  (14)  shows  that  the  computation  of 
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the  gradient  g(x)  requires  n additional  function  evaluations. 
In  turn,  if  f (x)  and  g(x)  are  known,  Eq.  (15)  shows  that  the 
computation  of  the  Hessian  H(x)  requires  m additional  function 
evaluations.  Clearly,  Eqs.  (14)  and  (15)  illustrate  the 
validity  of  Assumptions  (Al)  and  (A2).  Then,  providing  one 
can  ascertain  that  Assumption  (A3)  is  true,  one  concludes  that 
the  parameter  (11)  is  a correct  indicator  of  computational 
speed. 

Example  4.2.  Suppose  that  the  function  f depends  on  x, 
not  directly,  but  indirectly  through  some  variable  y, 
which  is  a function  of  x defined  by  means  of  some  definite 
integral.  For  simplicity,  assume  that  x and  y are  both 
n-vectors.  Then,  the  situation  is  as  follows: 

f = f (y)  , (16) 


where 


y 


,1 

<P  (x,  t)  dt 

0 


(17) 


Denote  by 


F (x)  = f (y  (x)  ) 


(18) 


k 
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the  function  obtained  by  combining  (16) -(17)  and  eliminating 
y.  Observe  that 


F=yf,  F =yf+yfy 
x xy  xx  -'xxy  x yy  x ' 


(19) 


Next,  assume  that  some  particular  algorithm  is  employed 

in  order  to  find  the  minimum  of  the  function  F (x) , utilizing 

the  gradient  (19-1)  and  perhaps  the  Hessian  (19-2).  Since 

the  computation  of  F,  F , F requires  previous  numerical 

X XX 

integrations,  defined  through  (17)  and  (20),  this  example 
illustrates  a situation  where  the  validity  of  Assumption 
(A3)  is  plausible.  Then,  providing  one  can  ascertain  that 
Assumptions  (Al)  and  (A2)  are  true,  one  concludes  that  the 
parameter  (11)  is  a correct  indicator  of  computational 
speed. 


17 


NPS-3 


5 . Cases  Where  the  Assumptions  Are  Not  Satisfied 

In  this  section,  we  present  some  cases  where  the  equivalent 
number  of  function  evaluations  (11)  cannot  be  regarded  to  be 
a correct  indicator  of  computational  speed. 

Example  5.1.  Suppose  that  the  function  (1)  is  such  that 
analytical  expressions  for  the  components  of  the  gradient  and 
the  elements  of  the  Hessian  matrix  are  available.  In  par- 
ticular, assume  that  the  function  (1)  has  the  quadratic  form 

f (x)  = a + bTx  + j xTcx  , (21) 

with  the  implication  that 

g(x)  = b + cx  , H(x)  = c , (22) 

In  (21)- (22),  a,b,c  are  constants  having  appropriate  dimen- 
sions. 

Next,  consider  the  operational  count  associated  with 

function,  gradient,  and  Hessian  evaluations.  Observe  that 

2 

the  computation  of  the  function  requires  (n+1)  multiplications 

and  n(n+l)  sums  and  that  the  computation  of  the  gradient 
2 2 

requires  n multiplications  and  n sums.  If  we  neglect  the 
addition  times  by  comparison  with  the  multiplication  times, 
we  arrive  at  the  following  conclusions  for  the  ratios  of  the 
basic  times  tq'ti/t2  : 


V_ 


0 


(23) 


Vt0  ~ n2/(n+l)2  , t2//t0  = 

Hence,  for  n relatively  large,  we  have 

Tl/T0  = 1 ' T2/T0  = ° • (24) 

Therefore,  it  appears  that,  from  the  point  of  view  of  the 
CPU  time,  one  gradient  evaluation  is  equivalent  to  one  (not  n) 
function  evaluation,  and  one  Hessian  evaluation  is  equivalent 
to  zero  (not  m)  function  evaluations.  As  a conclusion,  it 
appears  that  Assumptions  (Al)  and  (A2)  are  not  justified  for 
the  quadratic  function  (21). 

Remark.  The  result  (24-1)  represents  a worst-case  condition, 
because  function  evaluation  and  gradient  evaluation  have  been 
regarded  to  be  separate  operations.  Had  one  accounted  for 
the  commonality  of  the  product  cx  to  both  f (x)  and  g(x),  then 
Eqs.  (24)  would  have  been  modified  as  follows: 

Ti/t0  = 0,  t2/t0  = 0 , (25) 

thereby  invalidating  Assumptions  (Al)  and  (A2)  to  an  even 
larger  degree. 

Example  5.2.  Suppose  that  the  function  (1)  is  such  that 
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the  Hessian  matrix  has  a banded  structure.  For  example, 
consider  the  following  generalized  Rosenbrock  function 
(see  Oren,  Ref.  6) : 


n “ 1 rt  n-1  n ry 

f (x ) = I (x.-l)  + 100  Z ( x . - x.  ) , (26) 

i=l  1 i=l  1 


and  observe  that  the  associated  Hessian  matrix  is  tri- 
diagonal. Therefore,  only  those  elements  H^  which  are 
located  on  the  principal  diagonal  and  on  a contiguous  sub- 
diagonal need  to  be  computed,  since  all  of  the  remaining 
elements  vanish.  The  number  of  nonzero  elements  of  the 
Hessian  matrix  that  need  to  be  computed  is 


m = 2n-l  , 


(27) 


instead  of  m.  Table  1 shows  the  values  of  m,  m,  and  m/n  for 
n ranging  between  5 and  30.  Note  that  the  ratio  m/m  decreases 
as  n increases  and  becomes  of  order  1/10  for  n=30.  Therefore, 
even  if  finite-difference  methods  were  to  be  employed  in  the 
computation  of  the  Hessian  matrix,  it  appears  that  Assumption 
(A2)  is  not  justified  for  the  generalized  Rosenbrock  function 
(26)  . 
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Table  1.  Generalized  Rosenbrock  function. 


n 

in 

m 

m/m 

5 

/ 

9 

15 

0.600 

10 

19 

55 

0.345 

15 

29 

120 

0.242 

20 

39 

210 

0.186 

25 

49 

325 

0.151 

30 


59 


465 


0.127 
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6 . New  Definition  of  Equivalent  Number  of  Function  Evaluations 

Prom  the  examples  of  the  previous  sections,  it  appears 
that  there  are  cases  where  Eq.  (11)  cannot  be  regarded  as 
representative  of  the  computational  speed  of  an  algorithm, 
in  that  it  might  overestimate  the  importance  of  the  gradient 
contribution  and  the  Hessian  contribution  to  the  equivalent 
number  of  function  evaluations.  In  addition,  Eq.  (11)  dis- 
regards the  contribution  due  to  algorithmic  operations. 

In  an  attempt  to  correct  the  above  situation,  we  supply 
here  a new  definition  of  equivalent  number  of  function 
evaluations.  Specifically,  we  introduce  a more  general 
parameter  Ng  , which  is  constructed  so  as  to  reflect  ac- 
curately the  computational  effort  associated  with  function 
evaluations  and  algorithmic  operations. 

While  we  retain  Eqs.  (5 )-(8),  we  replace  the  assumptions 
expressed  by  Eqs.  (9)  with  the  following  definitions: 


(Bl) 

T1  = C1T0  ' 

(28-1) 

(B2 ) 

t2  = C2t0  ' 

(28-2) 

(B3) 

T = C,T 
a 3 e 

(28-3) 

Here,  C^,  C2'  C3  are  coefficients  to  be  determined  experimentally 
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or  theoretically  through  an  operational  count.  With  this 
understanding,  Eq.  (8)  can  be  rewritten  as 


T 


x.N 
0 e 


(29) 


where 


Ne  = (1  + C3)  (NQ  + C1N1  + C2N2)  . (30) 

This  expression  constitutes  a new  definition  of  equivalent 
number  of  functions  evaluations. 

Alternative  Definition.  If  one  introduces  the  new 
coefficients 


K1  = Cl/n  ' K2  = C2/m  ' K3  = 1+C3  ' 

Eq.  (30)  can  be  rewritten  in  the  alternative  form 

Ne  = K3 (Nq  + K1nN1  + K2mN2) 

Remark.  With  the  terminology  of  this  section, 
(Al),  (A2),  (A3)  can  now  be  restated  as  follows: 

(Al)  = n or  = 1 , 


(31) 


(32) 


Assumptions 


(33-1) 


¥ 


W 
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(A2 ) C2  = m or  = 1 , 


(33-2) 


(A3;  = 0 or  K,  = 1 


(33-3) 


Comment.  The  coefficients  and  C 2 measure  the  relative 
importance  of  the  computational  effort  associated  with  gradient 
evaluation  and  Hessian  evaluation  vis-a-vis  the  computational 
effort  associated  with  function  evaluation.  The  coefficient 
measures  the  relative  importance  of  the  computational 
effort  associated  with  algorithmic  operations  vis-a-vis  the 
computational  effort  associated  with  function,  gradient,  and 
Hessian  evaluations.  While  the  coefficients  and  C2  depend 
on  the  nature  of  the  function  f (x) , the  coefficient  depends 
also  on  the  structure  of  the  particular  algorithm  and  search 
technique  employed. 

The  above  coefficients  can  be  determined  through  either 
an  operational  count  or  through  computer  experimentation. 

By  determining  the  triplet  (C^,  C 2,  C^)  and  the  associated 
triplet  (K^ , , K^),  and  by  measuring  the  deviation  of  these 

coefficients  from  the  idealized  values  (33) , one  can  supply 
an  answer  to  the  basic  questions  formulated  in  this  paper, 
namely,  those  concerning  the  correctness  of  Assumptions  (Al) , 
(A2),  (A3). 
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7 . Numerical  Experiments 

In  the  following  sections,  we  describe  some  numerical 
experiments,  leading  to  the  computation  of  the  coefficients 
and  for  several  test  functions  and  minimization  al- 
gorithms . 

All  computations  were  performed  on  the  IBM  370/155  com- 
puter of  Rice  University.  FORTRAN  programming  was  employed 
in  conjunction  with  double-precision  arithmetic.  A FORTRAN 
G1  compiler  was  used.  A FORTRAN  TIME  subroutine  was  used 
in  order  to  determine  the  CPU  times.  Note  that  the  IBM 
370/155  computer  of  Rice  University  has  multi-programming 
and  time-sharing  capabilities. 

Test  Functions.  Five  test  functions  were  employed, 
namely,  the  Rosenbrock,  Wood,  Powell,  and  Miele  functions 
plus  the  generalized  Rosenbrock  function.  The  results  per- 
taining to  twelve  additional  test  functions  can  be  found  in 
Ref.  7. 

Example  7.1.  Rosenbrock  Function: 
f(x)  = (x^l)2  + 100  (x2  - x2)2  ; 

Example  7.2.  Wood  Function: 

f(x)  = (xL-l)2  + 100 (x2  - x,)2  + (x3-l)2  + 90 (x2  - x4)‘ 
+ 10.1  [ (x2-l ) 2 + (x4-l)2  J + 19.8(x2-1) (x4-1)  | 


(34  ) 


(35) 
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Example  7.3.  Powell  Function: 


f (x)  = (x1+  10x2)2  + 5(x3  - x4)2  + (x2  - 2x3)4 


+ 10 (x3  - x4) 


Example  7.4.  Miele  Function: 


4 6 

f(x)  = (exp  x^  - x2)  + 100  (x2  - x3) 


+ tan4  (x3  - x4)  + x®  + (x4  - l)2 


Example  7.5.  Generalized  Rosenbrock  Function: 


f(x)  = I (x^lp  + 100  Z (xf  - xi+1) 
i=l  i=l 


For  the  generalized  Rosenbrock  function,  values  of  n ranging 
between  5 and  30  were  considered. 

Minimization  Algorithms.  Three  unconstrained  minimiza- 
tion algorithms  were  employed:  The  Fletcher-Reeves  algorithm 

(FR) , the  Davidon-Fletcher-Powell  algorithm  (DFP ) , and  the 
quasilinearization  algorithm  (QL) . The  FR  and  DFP  algorithms 
exemplify  first-order  methods,  and  the  QL  algorithm  exem- 
plifies second-order  methods.  For  the  detailed  structure  of 
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these  algorithms,  see  for  example  Ref.  7. 

Fletcher-Reeves  and  Davidon-Fletcher-Powell  Algorithms. 
For  these  algorithms,  one-step  cubic  interpolation  was  em- 
ployed in  order  to  determine  the  optimum  stepsize  a.  In 
the  cubic  interpolation  process,  two  ordinates  and  two 
slopes  were  employed.  Therefore,  it  was  assumed  that  one 
iteration  of  each  of  these  algorithms  requires  two  function 
evaluations  and  two  gradient  evaluations.  For  the  details, 
see  Refs.  7-8. 

Quasilinearization  Algorithm.  For  this  algorithm,  the 
stepsize  a=l  was  employed.  Therefore,  it  was  assumed  that 
one  iteration  requires  one  function  evaluation,  one  gradient 
evaluation,  and  one  Hessian  evaluation.  The  solution  of  the 
linear  system  governing  the  components  of  the  search  direction 
was  done  by  Gaussian  elimination  without  pivoting.  For  the 
details,  see  Ref.  7. 

Coefficients  and  C These  coefficients  can  be 
computed  with  the  aid  of  Eqs.  (28-1)  and  (28-2)  as  follows: 

C1  = VT0  ' C2  = 1 2^T  0 ’ (39 

For  the  basic  times  Tq,  t^,  to  be  sufficiently  precise, 
it  is  necessary  that  the  computation  of  f(x),  g(x),  H(x)  be 
repeated  a large  number  of  times  at  the  same  nominal  point  x 
(for  instance,  1000  times) . 
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With  the  above  considerations  in  mind,  the  basic  times 
t , t ^ were  determined  as  follows:  (i)  by  evaluating 
separately  the  function,  the  gradient,  and  the  Hessian 
1000  times;  (ii)  by  monitoring  the  associated  CPU  times, 
which  include  the  do-loop  time;  (iii)  by  determining 
separately  the  do-loop  time;  and  (iv)  by  subtracting  the 
do-loop  time  from  the  experimentally  determined  CPU  times. 

With  the  basic  times  Tq,  t^,  T2  known,  the  coefficients 
C^  and  C2  can  be  computed  with  (39).  Then,  the  coefficients 
and  can  be  determined  with  (31-1)  and  (31-2). 

Coefficient  C^.  This  coefficient  can  be  computed  with 
the  aid  of  Eqs.  (5),  (6),  and  (28-3)  as  follows: 

C3  = (T-T0-T1-T2)/(T0+T1+T2)  • (40) 

Next,  we  invoke  Eqs.  (7)  and  the  definition 

T = tN  , (41) 


where  N denotes  the  number  of  iterations  and  t denotes  the 
time  required  to  perform  one  iteration  (this  includes  both 
the  algorithmic  time  and  the  function  evaluation  time) . In 
the  light  of  (7)  and  (41),  Eq.  (40)  can  be  rewritten  as 


C3  (tN-t0N0-t1N1-t2N2) / (t0N0+t1N1+t2N2) 


(42) 
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Examination  of  one  iteration  of  the  algorithms  under 
consideration  shows  that,  for  the  search  conditions  assumed, 


Nq/N  = 2 , Nj/N  = 2 , N2/N  = 0 (43) 

for  the  FR  and  DFP  algorithms  and  that 

Nq/N  = 1 , Nj/N  = 1 , N2/N  = 1 (44) 

for  the  QL  algorithm.  As  a consequence,  Eq . (42)  simplifies 

to 

C3  = (x  - 2tq  - 2t1)/(2tq  + 2t1)  (45) 

for  the  FR  and  DFP  algorithms,  and  to 

C3  = (T  “T0  “ Tl“  T2)/(T0  + T1  + x2)  (46) 

for  the  QL  algorithm.  For  the  basic  time  t to  be  sufficiently 
precise,  it  is  necessary  that  one  iteration  of  each  of  the 
algorithms  under  consideration  be  repeated  a large  number  of 
times  at  the  same  nominal  point  x (for  instance,  1000  times). 

With  the  above  considerations  in  mind,  the  basic  time  t 
was  determined  as  follows:  (i)  by  executing  one  iteration  of 
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each  of  the  algorithms  under  consideration  1000  times;  (ii) 
> by  monitoring  the  associated  CPU  times,  which  include  the 

do-loop  time;  (iii)  by  determining  separately  the  do-loop 
time;  and  (iv)  by  subtracting  the  do-loop  time  from  the  ex- 
perimentally determined  CPU  times. 

With  the  basic  time  t known,  and  with  Tq,  t^,  also 
known,  the  coefficient  C^  can  be  computed  with  (45)  for  the 
FR  and  DFP  algorithms  and  with  (46)  for  the  QL  algorithm. 
Then,  the  coefficient  can  be  determined  with  (31-3) . 
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8 . Numerical  Results 

Careful  numerical  experiments  were  performed  along  the 
lines  outlined  in  Section  7,  and  the  results  are  given  in 
Tables  2-5.  Tables  2-3  give  the  coefficients  C^,  C 2 and  K^, 

K2  for  the  test  functions  (34) -(38).  Tables  4-5  give  the 
coefficient  for  the  test  functions  (34) -(38)  and  the  FR, 
DFP,  and  QL  algorithms.  The  coefficient  is  not  given, 
since  it  can  be  computed  with  the  simple  relation  (31-3). 

Coefficient  C^.  Table  2 shows  that  the  coefficient  C^ 
for  the  Rosenbrock  function  has  the  value  1.04  (instead  of 
n=2) , and  the  coefficient  C^  for  the  Powell  function  has  the 
value  1.21  (instead  of  n=4).  Thereby,  use  of  the  standard 
definition  (11)  of  equivalent  number  of  function  evaluations 
overestimates  the  effort  associated  with  gradient  computation 
by  a factor  of  2 for  the  Rosenbrock  function  and  by  a factor 
of  3 for  the  Powell  function. 

Table  3 refers  to  the  generalized  Rosenbrock  function 
and  shows  that,  for  the  case  n=5,  the  coefficient  C^  has  the 
value  1.17  (instead  of  5);  for  the  case  n=30,  the  coefficient 
C^  has  the  value  1.28  (instead  of  30).  Thereby,  use  of  the 
standard  definition  (11)  of  equivalent  number  of  function 
evaluations  overestimates  the  effort  associated  with  gradient 
computation  by  a factor  of  4 for  n=5  variables  and  by  a factor 
of  24  for  n=30  variables. 


) 

1 
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In  conclusion,  under  the  hypothesis  that  analytical 
expressions  for  the  components  of  the  gradient  are  available, 
and  for  Examples  7.1  through  7.5,  it  does  not  appear  that 
Assumption  (Al)  is  satisfied. 

Coefficient  C 2.  Table  2 shows  that  the  coefficient  C2 
for  the  Rosenbrock  function  has  the  value  1.01  (instead  of 
m=3) , and  the  coefficient  C2  for  the  Powell  function  has  the 
value  0.98  (instead  of  m=10) . Thereby,  use  of  the  standard 
definition  (11)  of  equivalent  number  of  function  evaluations 
overestimates  the  effort  associated  with  Hessian  computation 
by  a factor  of  3 for  the  Rosenbrock  function  and  by  a factor 
of  10  for  the  Powell  function. 

Table  3 refers  to  the  generalized  Rosenbrock  function 
and  shows  that,  for  the  case  n=5,  the  coefficient  C 2 has  the 
value  1.42  (instead  of  m=15) ; for  the  case  n=30,  the  co- 
efficient C2  has  the  value  3.56  (instead  of  m = 465). 

Thereby,  use  of  the  standard  definition  (11)  of  equivalent 
number  of  function  evaluations  overestimates  the  effort 
associated  with  Hessian  computation  by  a factor  of  11  for 
n=5  variables  and  by  a factor  of  130  for  n=30  variables. 

In  conclusion,  under  the  hypothesis  that  analytical 
expressions  for  the  elements  of  the  Hessian  matrix  are 
available,  and  for  Examples  7.1  through  7.5,  it  does  not 
appear  that  Assumption  (A2)  is  satisfied. 


m>  - 
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Coefficient  C.^.  Inspection  of  Tables  4-5  shows  that  the 
coefficient  is  never  negligible  by  comparison  with  1, 
meaning  that  the  algorithmic  time  is  never  negligible  by 
comparison  with  the  function  evaluation  time.  Indeed,  in 
many  cases,  can  be  larger  than  1,  meaning  that  the  al- 
gorithmic time  can  be  larger  than  the  function  evaluation 
time. 

For  example,  consider  the  Wood  function.  Table  4 shows 
that  the  coefficient  has  the  value  1.14  for  the  FR  al- 
gorithm, 3.18  for  the  DFP  algorithm,  and  2.09  for  the  QL 
algorithm. 

As  another  example,  consider  the  generalized  Rosenbrock 
function.  Table  5 shows  that,  for  values  of  n ranging  be- 
tween 5 and  30,  the  coefficient  ranges  between  0.56  and 
0.73  for  the  FR  algorithm,  between  2.87  and  10.67  for  the 
DFP  algorithm,  and  between  1.94  and  22.05  for  the  QL  al- 
gorithm. 

In  conclusion, under  the  hypothesis  that  analytical  ex- 
pressions for  the  components  of  the  gradient  and  the  ele- 
ments of  the  Hessian  matrix  are  available,  and  for  Examples 
7.1  through  7.5,  it  does  not  appear  that  Assumption  (A3)  is 
satisfied  for  the  FR,  DFP,  and  QL  algorithms. 

Comments . With  reference  to  the  previous  results,  the 
following  comments  are  pertinent. 
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(i)  In  devising  the  subroutines  necessary  to  the  compu- 
tation of  the  functions  (34)- (38)  and  their  first  and  second 
derivatives,  an  effort  was  made  to  render  the  CPU  time  as 
small  as  possible.  The  commonality  existing  between  the  com- 
ponents of  the  gradient  was  exploited.  The  commonality 
existing  between  the  elements  of  the  Hessian  matrix  being 
computed  was  also  exploited.  Finally,  the  fact  that  the 
Hessian  matrix  might  have  many  elements  which  vanish  was  taken 
into  consideration. 

(ii)  For  the  subroutines  mentioned  in  (i) , the  computa- 

tion of  the  function,  the  computation  of  the  gradient,  and  the 
computation  of  the  Hessian  matrix  were  conceived  to  be  sepa- 
rate operations.  No  attempt  was  made  to  exploit  the  common- 
ality between  function  and  gradient  nor  the  commonality  between 
function,  gradient,  and  Hessian  matrix.  In  other  words,  a 
worst-case  situation  was  postulated,  and  conservative  esti- 
mates of  were  arrived  at.^  Had  the  above  commonality 

been  exploited,  the  values  of  C^,  C ^ would  have  been  smaller 
than  those  given  in  Tables  2-3  and  the  values  of  would 

have  been  larger  than  those  given  in  Tables  4-5. 


I 


Here,  the  adjective  conservative  is  employed  in  the  sense 
of  "favorable"  to  the  old  definition  (11)  of  equivalent 
number  of  function  evaluations. 


* 
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(iii)  For  the  generalized  Rosenbrock  function,  Table  5 
illustrates  in  a striking  way  the  effect  of  the  size  of  the 
problem  on  the  relative  importance  of  algorithmic  time  vis-a- 
vis  function  evaluation  time.  As  n increase,  tends  to  be 
constant  for  the  FR  algorithm;  it  increases  at  a linear  rate 
for  the  DFP  algorithm;  and  it  increases  at  a faster-than- 
linear  rate  for  the  QL  algorithm. 

The  explanation  for  this  result  is  quite  simple.  As  n 

increases,  the  function  evaluation  time  T increases  linearly 

e 

for  the  FR  and  DFP  algorithms  and  quadratically  for  the  QL 
algorithm.  On  the  other  hand,  as  n increases,  the  algorithmic 
time  T increases  linearly  for  the  FR  algorithm,  quadratically 

d 

for  the  DFP  algorithm,  and  cubically  for  the  QL  algorithm. 

(iv)  The  coefficients  and  depend  mostly  on  the 
nature  of  the  function  f(x).  On  the  other  hand,  the  coef- 
ficient for  the  FR  and  DFP  algorithms  depends  also  on  the 
search  technique  employed  to  determine  the  stepsize  a;  should 

a different  search  technique  be  employed,  the  value  of  would 
change.  By  the  same  token,  the  coefficient  for  QL  algorithm 
depends  on  the  subroutine  used  to  solve  the  linear  system 
governing  the  components  of  the  search  direction;  should  a 
different  subroutine  be  employed,  the  value  of  would  change. 

(v)  In  this  paper,  the  coefficients  C^,  C 2,  C3  were 
determined  through  computer  experiments.  An  alternative  way 


is  to  determine  these  coefficients  through  an  operational 
count.  This  operational  count  requires  the  a priori  know- 
ledge of  the  times  necessary  to  perform  each  computer 
operation;  in  other  words,  it  requires  the  a priori  knowledge 
of  the  relative  weight  of  each  computer  operation.  An 
analysis  along  these  lines  has  been  carried  out  in  Ref.  7 
and  leads  to  results  qualitatively  consistent  with  those 
presented  here.  However,  this  analysis  is  omitted  for  the 
sake  of  brevity. 
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Table  2.  Results  for  the  Rosenbrock,  Wood,  Powell,  and 
Miele  functions. 


Example  n m 


3 

1.04 

1.01 

0.52 

0.34 

10 

1.02 

0.81 

0.26 

0.08 

10 

1.21 

0.98 

0.30 

0.10 

10 

1.29 

1.78 

0.32 

0.18 

Table  3.  Results  for  the  generalized  Rosenbrock  function. 
Example  n m 


. 5 

5 

15 

1.17 

1.42 

0.23 

0.09 

.5 

10 

55 

1.20 

1.82 

0.12 

0.03 

. 5 

15 

120 

1.23 

2.36 

0.08 

0.02 

.5 

20 

210 

1.21 

2.71 

0.06 

0.01 

.5 

25 

325 

1.28 

3.15 

0.05 

<0.01 

. 5 

30 

465 

1.28 

3.56 

0.04 

<0.01 
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Table  4.  Results  for  the  Rosenbrock,  Wood,  Powell, and  Miele 
functions . 


Algorithm 

FR 

DFP 

QL 

Example 

n 

m 

C3 

C3 

C3 

7.1 

2 

3 

1.65 

2.68 

1.17 

7.2 

4 

10 

1.14 

3.18 

2.09 

7.3 

4 

10 

1.28 

3.31 

2.21 

7.4 

4 

10 

0.43 

1.05 

0.55 

Table  5.  Results 

for  the 

generalized 

Rosenbrock 

function . 

Algorithm 

FR 

DFP 

QL 

Example 

n 

m 

C3 

C3 

C3 

7.5 

5 

15 

0.73 

2.87 

1.94 

7.5 

10 

55 

0.68 

4.24 

4.89 

7.5 

15 

120 

0.62 

5.70 

8.09 

7.5 

20 

210 

0.61 

7.46 

12.54 

7.5 

25 

325 

0.58 

9.01 

17.34 

7.5 

30 

465 

0.56 

10.67 

22.05 
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9 . Conclusions  and  Recommendations 

In  this  paper,  we  have  considered  the  comparative 
evaluation  of  algorithms  for  mathematical  programming 
problems  from  the  point  of  view  of  computational  speed.  We 
have  examined  critically  the  indirect  measurement  of  com- 
putational speed  through  the  equivalent  number  of  function 
evaluations  Ng,  and  we  have  found  that  the  Ng  concept, 
while  accurate  in  some  cases,  has  drawbacks  in  other  cases. 
Indeed,  it  might  lead  to  a distorted  view  of  the  relative 
importance  of  an  algorithm  with  respect  to  another. 

In  an  effort  to  correct  the  above  distortion,  we  have 
imbedded  the  parameter  into  a more  general  parameter 
Ng  , which  is  constructed  so  as  to  reflect  accurately  the 
computational  effort  associated  with  function  evaluations 
and  algorithmic  operations.  This  new  parameter  Ng  in- 
cludes coefficients  C^,  which  can  be  determined 

either  experimentally  or  through  an  operational  count.  When 
the  triplet  (C^,  C C^)  takes  on  the  values  (n,m,0),the 
new  parameter  Nq  reduces  to  the  old  parameter  Ng . 

We  have  determined  experimentally  the  coefficients 
C^,  C2 , for  five  test  f’  actions  and  three  minimization 
algorithms.  And  we  have  found  that  the  deviations  of  these 
coefficients  from  the  idealized  values  n,m,0  can  be  sub- 


stantial. 
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More  specifically,  the  experimental  values  found  for 
are  2 to  24  times  smaller  than  the  idealized  value  n. 

The  experimental  values  found  for  C ^ are  3 to  130  times 
smaller  than  the  idealized  value  m.  And  the  experimental 
values  found  for  are  never  negligible  with  respect  to 
1;  they  are  of  order  1 for  the  FR  algorithm,  and  of  order  1 to 
10  for  the  DFP  and  QL  algorithms. 

Obviously,  the  experimental  values  of  C^,  C are 
subject  to  errors  due,  among  other  things,  to  the  multi- 
programming and  time-sharing  capabilites  of  the  IBM  370/155 
computer  of  Rice  University.  However,  the  basic  fact  re- 
mains that  the  deviations  detected  for  C^,  , C-,  from  the 

idealized  values  n,m,0  are  so  large  that  the  use  of  the  Ng 
concept  is  open  to  serious  question.  By  the  way,  conclusions 
along  these  general  lines  have  also  been  obtained  by  Hillstrom 
in  Refs.  9-10. 

From  the  analyses  performed  and  the  results  obtained,  it 
is  inferred  that,  due  to  the  weakness  of  the  Ng  concept,  the 
use  of  the  Nq  concept  is  advisable  as  a means  for  comparing 
different  algorithms  from  the  point  of  view  of  computational 
speed.  In  effect,  this  is  the  same  as  stating  that,  in  spite 
of  its  obvious  shortcomings,  the  direct  measurement  of  the 
CPU  time  is  still  the  most  reliable  way  to  compare  different 
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minimization  algorithms. 

However,  for  the  direct  measurement  of  the  CPU  time  to 
be  really  meaningful,  provisions  similar  to  those  implemented 
in  Refs.  1-2  should  be  employed.  That  is,  it  is  necessary 
that  the  comparison  of  different  algorithms  be  done  on  a 
single  computer,  with  the  same  programming  language,  with 
the  same  compiler,  with  the  same  subroutines,  under  similar 
workload  conditions  of  the  computer,  and  by  the  same  pro- 
grammer. 

In  closing,  these  authors  stress  that  the  conclusions  of 
this  paper  should  not  be  interpreted  as  an  invitation  to 
other  authors  to  disregard  reporting  on  number 
of  iterations  N,  number  of  function  evaluations  Ng,  number 
of  gradient  evaluations  N ^ , and  number  of  Hessian  evaluations 
. By  all  means,  these  are  useful  quantitites,  which  should 
be  reported  because  their  knowledge  does  shed  some  light  on 
the  comparative  behavior  of  different  algorithms.  Neverthe- 
less, none  of  the  methods  for  the  indirect  measurement  of 
the  computational  speed  is  truly  satisfactory, and  these 
authors  feel  that  there  exist  no  reliable  alternative  to  the 
direct  measurement  of  the  CPU  time. 


* 
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